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By employing a nonlinear quantum kicked rotor model, we investigate the transport of energy in 
multidimensional quantum chaos. This problem has profound implications in many fields of science 
ranging from Anderson localization to time reversal of classical and quantum waves. We begin our analysis 
with a series of parallel numerical simulations, whose results show an unexpected and anomalous behavior. 
We tackle the problem by a fully analytical approach characterized by Lie groups and solitons theory, 
demonstrating the existence of a universal, nonlinearly-enhanced diffusion of the energy in the system, 
which is entirely sustained by soliton waves. Numerical simulations, performed with different models, show 
a perfect agreement with universal predictions. A realistic experiment is discussed in two dimensional 
dipolar Bose-Einstein-Condensates (BEC). Besides the obvious implications at the fundamental level, our 
results show that solitons can form the building block for the realization of new systems for the enhanced 
transport of matter. 

Anderson localization is a fundamental concept that, originally introduced in solid-state physics to describe 
conduction-insulator transitions in disordered crystals, has permeated several research areas and has 
become the subject of great research interest 1 " 11 . Theories and subsequent experiments demonstrated that 
disorder favors the formation of spatially localized states, which sustain diffusion breakdown and exponentially 
attenuated transmission in random media 1 . Although many properties of wave localization are now well under- 
stood, several fundamental questions remains. Perhaps one of the most intriguing problem is related to the 
transport of energy. Intuitively, one can expect that disorder -by favoring exponentially localized stated- arrests 
in general any propagation inside a noncrystalline medium. However, the interplay between localization and 
disorder is nontrivial 5,12,13 and under specific conditions randomness can significantly enhance energy transport. 
In linear regime, in particular, it has been observed that quasi-crystals with multifractal eigenstates and/or 
material systems with temporal fluctuations of the potential (or refractive index), lead to anomalous diffusion 
in the phase space 14 " 18 . This originates counterintuitive dynamics including ultralow conductivities 14 , as well as 
the formation of mobility edges even in one dimensional systems 17 . All these studies focused on linear materials 
and they did not investigate the role of nonlinearity for further controlling the transport of energy in many 
dimensions. 

In the context of quantum localization, the problem of energy transport has stirred a conspicuous interest as 
well. In this area, quantum-classical correspondences mediated by Anderson localization possess many implica- 
tions in the irreversible behavior of time reversible systems, which are at the basis of a long standing physical 
debate -i.e., the Loschmidt paradox 19 . Started a few centuries ago as a controversy between Boltzmann and 
Loschmidt, this famous paradox is still the matter of intense research in the scientific community 20 " 26 . The 
Loschmidt paradox deals with the origin of the irreversible behavior of time reversible systems that, according 
to the law of classical mechanics, should not manifest irreversible entropy growth as conversely predicted by the 
second law of thermodynamics. Time reversibility, in fact, guarantees that for every orbit that leads to an entropy 
increase there exists -with the same probability- a time-reversed path that generates the same entropy change but 
with an opposite sign 19 . Recent experiments might suggest that a possible solution of this paradox can be 
formulated in terms of deterministic chaos 23 . According to this interpretation, time reversibility is possible only 
at the quantum level, where Anderson localization breaks diffusive transport and suppresses the mixing ability of 
chaos, as discussed in one dimension by a series of papers 21,22,27 . The possibility to exploit quantum time reversal is 
at the basis of several nontrivial dynamics including the quantum-echo effect 21,22,26 . However, when more dimen- 
sions are considered, numerical simulations predict that ergodicity is fully restored and diffusive transport settles 
is again, thus re-establishing the classical features of chaos and preventing quantum time reversal and its 
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associated dynamics 22 . Nevertheless, theoretical work reported to 
date considered only noninteracting systems, characterizesd by lin- 
ear equations of motion. The Loschmidt paradox, conversely, 
involved the use of interacting atoms, whose interplay in the mean 
field regime is accounted by short and/or long ranged nonlinear 
responses 28 " 30 . As pointed out in the literature 21 , atoms interactions 
are of crucial importance in quantum localization and diffusion. A 
key question in this problem lies in understanding how nonlinearity 
affects the transport of energy in many dimensions. In one dimen- 
sional quantum chaotic systems, pioneering numerical experiments 
of Casati et al. 27 reported that nonlinearity maintains dynamical 
localization effects. This conforms to the intuitive idea that a non- 
linear response, due to its localization properties, works together 
Anderson effects to suppress mixing dynamics. Recent theoretical 
work performed on d- dimensional disordered lattices show that at 
any finite nonlinearity there exist a finite probability for the obser- 
vation of Anderson localization effects 31 . In this scenario, nontrivial 
effects are expected to occur on the energy transport, due to the rich 
interplay between localization and nonlinearity, as well as by the 
additional degrees of freedom that can interact in the dynamics. 

In this Article, we theoretically investigate this problem by 
employing both numerical simulations and analytic techniques. To 
pursue a general theory, we here consider the following two dimen- 
sional model: 



j^(r'-r)|iA(r')| 2 + ^rW 



= 0, 



(1) 



with r = (jc, y\ V 2 = d 2 /dx 2 + d 2 /dy\ 5 T = E„ S(t - nT) a periodic 
delta- function of period T, R a general nonlinear response and 
U(x,y) =y(cos x+ cos y)+e cos (*+;/) a two dimensional peri- 
odic potential with strength defined by e and y. Equation (1) defines 
a two dimensional, nonlinear quantum kicked rotor: for R = 0 it 
reduces to the linear quantum kicked rotator 22 while for U = 0 it 
corresponds to the 2D nonlinear Schrodinger equation (NLS), which 
represents a universal model of nonlinear waves in dispersive 
media 30 . In one dimension, conversely, Eq. (1) generalizes the non- 
linear model investigated in 27 with classical chaos parameter K = 
2yT. Despite its deterministic nature, Eq. (1) can be precisely mapped 
to the Anderson model with a random potential 11 ' 32,33 , and therefore 
furnishes a fundamental model for studying energy transport and 
dynamical localization 34 ' 35 in random systems. 

This article is organized as follows. We begin by a series of experi- 
ments studying the diffusion of energy by integrating Eq. (1) in time 
with a parallel algorithm. We interpreted the results by an analytic 
theory based on Lie symmetry groups and soliton waves, which 



predicts a universal enhanced diffusion process, which is entirely 
sustained by nonlinearity. Analytical results are validated against 
numerical simulations, snowing a perfect agreement with our pre- 
dictions. We finally discuss the realization of a possible experiments 
in a multidimensional dipolar BEC, showing the feasibility of our 
ideas in a realistic setting. 

Results 

A first numerical experiment. We begin our analysis by calculating 



the momentum diffusion (P) = 



\jj \ versus time, with p = V/i 



the momentum operator and = J drj\\j/\ 2 the quantum 

average. In order to investigate a general system, we considered a 
nonloncal diffusive nonlinear response n = j drR(r' — r)\i//(r')\ 2 
following from: 

(l-a 2 V 2 )n = \4,\ 2 , (2) 

with nonlocality controlled by a. When a = 0, the system response is 
local with n = \i//\ 2 . For a 0, conversely, the system nonlinearity 

becomes long ranged with kernel given by£( r ) = — K 0 - , being 

K 0 the modified Bessel function of second kind. Diffusive nonli- 
nearities are particularly interesting in the context of nonlinear 
optics, as they can be easily accessed in liquids, as well as in Bose- 
Einstein Condensates (BEC), where they generalize previously 
investigated models 36 " 38 . Parallel numerical simulations are per- 
formed by a direct solution of (1) with an unconditionally stable 
algorithm. In order for the field \jj to explore the periodic potential 
(7, we here consider wave packets whose spatial extension Ar <C 2n. 
Figure 1 summarizes our results obtained for a = 0.2, by launching at 
the input a gaussian beam i/f = Ae _ *V ca o with waist co 0 = 0.3 and 
amplitude A = 4 (Fig. la). The stochastic parameter K has been set to 
K= 1.8 > K*, above the stochastization threshold iC* ~ 0.97 where 
the linear classical uncoupled rotor exhibits diffusive transport in 
momentum space 22 . For comparison, we also calculated the linear 
dynamics resulting from R = 0 (Fig. lb dotted line). As seen from 
Fig. lb, the 2D nonlinear rotor behaves dramatically different with 
respect to its linear counterpart, demonstrating the strong role 
played by nonlinearity in the process. In particular, the linear 
system exhibits Anderson localization and diffusion suppression 
for e = 0 (uncoupled condition), while for growing e it shows a 
monotonically increasing sub-diffusion (Fig. lb). In the nonlinear 
regime, conversely, Anderson localization is suppressed even for 
e = 0, and the dynamics shows an erratic, random-like behavior 
that does not manifest any simple monotonic increase for growing 
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values of e. These results are also significantly different from the 
nonlinear kicked rotor in one dimension 27 , where nonlinearity was 
observed to induce Anderson localization effects thus suppressing 
any diffusive (or sub diffusive) regime. 

Analytic theory and Universal diffusion scaling. To theoretically 
investigate the results of Fig. 1 and derive predictions of universal 
character, we begin by calculating the time evolution of the wave 
packet center of mass, thus generalizing the Ehrenfest theorem of 
classical quantum mechanics: 



dt 



= 2Re{(*A|i#}}, 



dt 



= MViW, (3) 



with F = U3 T + j dr'R(r - r') |i/Kr') | 2 . Equations (3) are valid for any 
dimensionality of the problem and for any nonlinearity R. In the 
following we assume a generic nonlinear response R that support 
at least a stable bound state of minimum energy, i.e., a nonlinear 
ground state solution. Numerical simulations performed in the 
previous section showed that the spatial field profile of the 
wavepacket, despite the chaotic motion, is not significantly altered 
in time (Fig. la,c). We can therefore describe the wavepacket 
dynamics in terms of a reduced set of coordinates modeling the 
nonlinear ground state of Eq. (1). To find the general form of the 
nonlinear ground state, we exploit Lie symmetry groups theory 39 . In 
particular, we start from the Lagrangian density £ of Eq. (1), written 
for U=0: 



~df 



dt 



J?(r-r')|<Kr')| 2 + |V^| 2 , (4) 



and identify its variational symmetries, which we express by the 
following basis of Lie generators: 



vi=V, 



v 3 -- 



d_ 

It' 



dx 



v 4 = ixjj 



(5) 



These generators are associated to translational, rotational and 
Gauge symmetries of (1). The nonlinear ground state of (5) 
represents an invariant solutions with respect to the global 
symmetry group generated by: 

v = r 0 -vi + v 2 + v 3 + /?v 4 . (6) 

We employed the method of characteristics 39 ' 40 to find the functional 
form of the general solution, which reads as follows: 



•A = <H|r-r 0 |y' 



(7) 



with <j> being a complex envelope. Equation (7) represents a soliton 
wave of the system. Closed form expression of </> for integrable 
responses R are found by the inverse scattering transform. In the 
case of Eq. (2), which possesses a nonintegrable response, we 



found approximate solutions by a variational analysis 36 ' 41 . In 
particular, we use the following Gaussian ansatz: 

V 71 a 

defined by the power P = (i// 1 1//) and waist a(t) . By substituting Eq. (8) 
in (4), after long but straightforward algebra we obtain the classical 
dynamics following from the Hamiltonian Tt\ 



+v, v=- 



2na 2 \Ao 2 



(9) 



with Z(x) = e~ x T(0,x), T(0,x) the upper incomplete gamma function 
and V the potential of the one dimensional motion of a. The potential 
V has a bell shape profile that possesses a unique absolute minimum 
V(a*) for every combination of P and a. The fixed point a(0) = a* 
corresponds to a soliton wave of the system, which propagates in a 
translational fashion with fixed waist a(i) = a*, while different initial 
values lead to a breather 42 characterized by a periodic oscillation of a 
in time. Figure 2a shows a typical V profile obtained for P = 100 and 
different nonlocality degrees a. As seen in the figure, the potential 
well becomes bigger and bigger for increasing non locality, 
explaining the dynamical robustness of the soliton dynamics 
during the chaotic motion observed in the parallel simulations 
performed in the previous section. 

To investigate the motion of the soliton ground state in the general 
case when U # 0 and R is arbitrary, we substitute Eq. (7) into (3) and 
perform an integration from nT to (n + 1)T. After some algebra, we 
end with the following system: 

Vn+l=Vn-[ K -g + 2eT Sm(Xo+yoH, 

q«+i = q»+p»+i> 

with g = [sin x 0 , smy 0 ], u = [1, 1], classical position q„ = q(nT) and 
momentum p n = p(nT) defined from q = (i^|r 0 |i^)/(i^|i/^) and 
p = Re{(i/f|p|i/f)}-2T/(i/f|i/f), respectively. In the derivation of (10) 
we assumed the general condition of (i^|(r — r 0 ) 2 |^) <ti 2n, in agree- 
ment with our introductory premises. Equations (10) represent a two 
dimensional standard map: for K > K*, above the stochastization 
threshold of the single uncoupled rotor, Eq. (10) is hyperchaotic and 
each dimension acts as an external noise source to the other, increas- 
ing the mixing of the overall system 43 . To highlight such a dynamics, 
we plot in Fig. 2b the positive Lyapunov exponent X calculated for 
Eqs. (10). As seen, even when only a single rotor overcomes its 
stochastization threshold, it provides a noise source to the other 
giving rise to two positive Lyapunov exponents. The largest 
Lyapunov exponent grows linearly with e (Fig. 2b). The hyperchaotic 
nature of Eq. (10) is expect to strongly affect the momentum dif- 
fusion in the phase space. We investigate the latter by generalizing 
the approach of Rechester and White developed for turbulent flows 44 . 
The classical diffusion D of the map can be expressed as follows: 



500 



> 0 



-500 




Figure 2 | (a) soliton potential V calculated for P = 100 and for different non locality degrees a; (b) positive Lyapunov exponent X of the chaotic map (10) 
versus coupling e, calculated for K = 5 and e = 0.8. 
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D= 



Ap 2 



^ p2tz p+go 

n Jo J _ x 



Ap 2 

2 ' 



(11) 



being Ap = p„ — p 0 and V the probability distribution of position p„ 
and momentum q n at time n (measured in kick units). To calculate V, 
we begin by observing that the the evolution of the map ( 10) at the n- 
th kick can be written as follows: 



Pn=Po + T "-l> <ln=<ln-l+Po+ T "-l> 



(12) 



expressed as a function of the initial momentum p 0 , with 
T h= Y?i=\ K '8j + 2eT sin(xy+7 ; -)-u and g ; = [sin x p sin y,]. The 
evolution of the probability density V at time n is then given by: 

+ go + go p+oo p2ti 

^(q„,p„,«)= ^ • • XM dp ° rfq ° x 

U,_ = _ rn J 1] = _ m J-GO Jo 



h n — - co h\ = — go 1 

P(qo,P 0 ,0)<5(p M -Po-T„-i) 



(13) 



being T^qo, po >0) the initial density. The presence of the additional 
summations over h originates from the periodicity of the position q 
in the phase space. By considering a general uniform initial distri- 
bution with all the particles possessing nonzero momentum p, i.e., 

^(qo'Po' 0 ) = — — — ' after substituting into (13) and (15), we 



obtain: 



271 



i ~r go + go c zn da- / x 

D=limi V ••• Y n ^T 2 „ ^(q,-q,- 1 -Po-T,- I ). 

"-»»,„e ro fc^^-Oo 4k - 1 (14) 

The leading order in the integral is obtained by neglecting any folding 
in the position space, i.e., by assuming hj = 0 (for hj 0, conversely, 
we get asymptotically small corrections expressed in terms of Bessel 
functions that we do not reported here due to their cumbersome 
expressions). The leading diffusion reads as follows: 

K 2 + Z 2 



D = 



15 



being ^ = 2eT. Equation (15) is to be considered of universal applic- 
ability, i.e., for any system dimension and for arbitrary nonlinear 
responses, as it has been derived under the general equations (3), 
(10)-(15). 

Discussion 

Equation (15) allows to fully interpret the nonlinear dynamics of Eq. 
(1). In particular, the quantum average (P) results from an hyperch- 
aotic system described by a two dimensional dimensional standard 



map, and each realization manifests itself as a random walk in Fig. lb. 
The map diffusion rate is identical to the momentum diffusion of the 
classical linear rotor 22 , hence, an additional average (in time or over 
an ensemble of input conditions) re-establishes a perfect classical 
correspondence for every coupling e > 0. It is worthwhile observing 
that the classical correspondence in the multidimensional linear 
quantum rotor is manifested only for very high coupling e, and in 
general the quantum diffusion (P) follows a fractional behavior with 
(P) oc t p<l (see e.g. 22 , or Fig. lb dashed lines). As a result, the linear 
quantum rotor sub-diffuses at a slower rate than its classical counter- 
part. Conversely, Eq. (15) predicts a perfect classical correspondence 
for every coupling e, which is re-established thanks to nonlinear 
effects. In order to verify Eq. (15), and to demonstrate such a non- 
linearly- enhanced transport dynamics, we performed extensive 
numerical simulations from Eq. (1) and calculated the average dif- 
fusion through a quantum average followed by an average over dif- 
ferent input conditions 



(p>= 



(16) 



Figure 3 summarizes our results obtained for K = 5, a = 0.2 and by 
considering an initial wave packet composed by a Gaussian beam 
with waist co 0 = 0.3 and amplitude A = 4. In complete agreement 
with Eqs. (10) -(15), we observe a diffusive behavior (P) oc t for every 
£>0 (Fig. 3a solid lines), whose rate is exactly matching our theor- 
etical prediction based on Eq. (15) [Fig. 3a dashed lines]. Figure 3b 
compares linear and nonlinear evolution of Eq. (1), the former 
obtained for R = 0. As seen in the figure, the nonlinear enhancement 
in the diffusion due to the restoration of classical effect is significant, 
and it increases with the coupling e. This can be intuitively expected 
due to the fact that an increasing coupling leads to an increasing 
mixing in the system and to stronger diffusion effects. The relative 
diffusion variation AD= DM0 "^ eflr ~ Dfaear , normalized respect to the 

^nonlinear 

nonlinear dynamics, is approximatively 50% at e = 1 .2, which means 
that nonlinearity is significantly faster (by a factor of two) in trans- 
porting energy. 

In order to discuss a possible experimental realization, and further 
verify the universality of our predictions, we considered the case of a 
multidimensional dipolar BEC 41 . This system attracted a conspicuous 
interest in the scientific community due to its important implications 
in many-body dynamics, quantum computing and nonlinear 
waves 28 ' 45 " 48 . Dipolar BEC are characterized by long range interactions, 
which support stable ground state solitons and high order azimu- 
tons 28,48 . The nonlinear dynamics of the wavefunction in a kicked 
optical lattice can be written in the following adimensional form: 

dk 



47I 2 



■h(k)h 



¥ = 0, (17) 




Figure 3 | (a) Momentum diffusion (P) versus time calculated from Eq. (1) (solid lines) and Eqs. (15) (dashed lines). Panel (b) compares the nonlinear 
diffusive dynamics based on Eq. (1) to the linear transport dynamics obtained in the case of R = 0. 
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Figure 4 | (a) Momentum diffusion (P) versus time calculated for a dipolar BEC in a kicked optical lattice (solid lines) and theoretical prediction based on 
the universal Eq. (15) (dashed lines). Panel (b) illustrate the spatial density probability distribution l^l 2 at time z = 50T. 



with dimensionless t = tco?, X = 



2mw z 



2mw z 



gPo 



where V = 



dX 2 + dY 2 ' 



P = 



h(k)=2-3nke k erfc(fc), 



is the Fourier transform of 



the BEC trapping frequency along z, I = A 



, g the coupling con- 



stant characterizing two -body contact interactions, g d the strength of 
the dipole-dipole interaction potential 28 and poo a reference BEC den- 
sity. Interaction terms g and gd can be tuned via Feshbach resonance, 
which allows to experimentally control the ratio g d lg that determines 
the properties of the nonlinear ground state of the condensate. In our 
numerical experiments, we considered g d lg = —0.5, as a feasible value 
for the generation of a stable two-dimensional soliton state 28 . 

Figure 4 summarizes our results for K = 2, T = 1 and different 
coupling strengths e. In perfect agreement with our theory, the aver- 
aged momentum diffuses as predicted by our universal formula, 
increasing its strength as the coupling e grows. The final density 
distribution, after t = 50 T, shows the stable propagation of the 
dipolar ground state that is evolving inside the system (Fig. 4b), 
proving the feasibility of our predictions in a realistic system. 

In conclusion, motivated by the large interest in the study of the 
energy transport in disordered system, we investigated the role of 
nonlinearity in affecting the dynamics of energy diffusion in random 
media. This problem has also profound implications in quantum 
chaos and time reversibility of classical systems, where nonlinearity 
emerges naturally when we consider interacting particles. We 
considered a universal model of wave propagation, namely the multi- 
dimensional quantum kicked rotor in presence of a generic atom- 
atom interactions. We began our analysis with a single numerical 
experiment, and studied the behavior of the energy diffusion in the 
phase space, finding unexpected results that are very different with 
respect to both linear multidimensional dynamics and nonlinear 
evolutions of one dimensional systems. We tackled the problem by 
employing a combination of soliton theory and Lie symmetry ana- 
lysis, finding a universal diffusion evolution of the energy that follows 
from the full restoration of classical effects sustained by soliton 
waves. The latter, in particular, breaks Anderson localization effects 
and diffuse energy with a larger rate with respect to linear systems of 
the same size. Numerical simulations performed on different models, 
including dipolar BEC, perfectly agree with our universal predic- 
tions. From a pure quantum perspective, our results demonstrate 
that atom-atom interactions inhibit quantum time reversal in many 
dimensions, due to the full recovery of classical chaotic mixing in the 
system. This generalizes the intuition of Adachi, Today and Ikeda 22 , 



who found that in linear regime multidimensional quantum time 
reversal is only conditionally possible, in the sense that an e threshold 
exists for the recovery of the initial wavepacket, while beyond a 
specific interaction e no reversal is possible. When nonlinearity is 
taken into account, conversely, solitons reestablish a fully mixing 
dynamics and for all values of e no time reversal can be observed. 

We can therefore conclude that nonlinear waves can favor the 
energy transport in a disordered medium, and significantly speed 
up the process of energy diffusion when compared to linear 
dynamics. Solitons can be therefore used for the development of 
new architectures for enhancing the transport of matter in disor- 
dered materials. 

Methods 

Numerical simulations of Eq. (1) have been realized by an homemade parallel code 
based on an unconditionally stable, second order time marching scheme. 
Parallelization is achieved by a two dimensional domain decomposition strategy, 
where each part of the computational domain is assigned to a different processor, with 
all communications written following the MPI standard. The numerical results pre- 
sented in this work have been realized by 200000 single cpu hours on 128 processors 
of our "reddragon" linux cluster. 
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